#include <math.h>
#include <stdio.h>
#ifndef FFT_H
#include "fft.h"
#endif // FFT_H
#ifndef FFTCONST_H
#include "fftconst.h"
#endif // FFTCONST_H
#include "fft.c"
#include "gettype.c"
#include "fft_intN.c"
#include "fftfixed.c"
#define N1 256
#define N 2048
#define PRINTSINE //定义了打印定点计算正弦数值
#ifdef PRINTSINE
#include "printfixedsin.c"
extern void printfixedsin(short n1,short n,short flag);
#endif
extern int fix_fft(myint x[], myint y[], short m, short inverse);
extern int fix_fftr(myint f[], short m, short inverse);

extern void fix_FFT(myint *x, myint *y, int n, int sign);
extern void FFT(myint *x, myint *y, int n, int sign);
extern void IFFT(myint *x, myint *y, int n, int sign);
extern void fft(int n, complex f[]);	//傅立叶变换 输出也存在数组f中
extern void ifft(int n, complex f[]);	// 傅里叶逆变换

int
fftsamples(myint x[],myint y[],int n1,
           int sample_pts_num, float Ampltitude, float frequence, float phi)
{
//N点FFT和IFFT,double x[],double y[]分别对应原序列的实部和虚部,调用FFT格式：FFT(x,y,N,1);
    /* int x[N] = {21949,20137,2875,-4142,4987,12483,7432,0,222,-1659,-13491,-24142,-15603,5995,14296,0, -14296 ,  -5995 ,   15603  ,  24142  ,  13491,1659  , 222  ,0  , -7432 ,  -12483 ,  -4987  ,  4142,-2875 ,  -20137  , -21949,  0};
       int y[N] = {0};

    myint x[N] = { 1, 5, 4, 8, 4, 1, 10, 14, 10, 11, 5, 7, 9, 1, 12, 3, 5, 14, 2, 2, 1, 8, 4, 5, 11, 10, 13, 2, 5, 9, 8, 6, 7, 1, 2, 3, 6, 5, 9, 8 };
    myint y[N] = { 0 };*/
    complex f[N];

    int i = 0;
    int n=1<<n1;

    for (i = 0; i < n; i++)
    {
        if (i < sample_pts_num)
        {
            x[i] =
                1 + sin(2 * PI * frequence * i / n +
                        phi) * Ampltitude +
                sin(2 * PI * 2 * frequence * i / n +
                    PI * 45 / 180) * 2 * Ampltitude;
        }
        else
        {
            x[i] = 0;
        }
        y[i] = 0;
    }

    return 0;
}

int printxy(FILE*fp,int n,myint x[],myint y[])
{
    int i;
        for (i = 0; i < n; i++)
    {
        /*printf("%6d %6d %6d\n", i, x[i], y[i]);
        fprintf(fp, "%6d %6d %6d\n", i, x[i], y[i]);*/
        printf("%6d %g %g\n", i, x[i], y[i]);
        fprintf(fp, "%6d %g %g\n", i, x[i], y[i]);
    }
    return 0;
}

int
callfft(myint x[],myint y[],unsigned char imethod,short inverse,unsigned char n1,unsigned char m,float Ampltitude, float frequence, float phi)
{

//N点FFT和IFFT,myint x[], y[]分别对应原序列的实部和虚部,调用FFT格式：FFT(x,y,N,1);
    myint z[N];//Z存放x,y排序后的索引，特征频率索引从大到小
    complex f[N];
    int i;

    unsigned int n=1<<n1;
    unsigned int sample_pts_num=1<<m;
    FILE *fp;
    fp = fopen("fft.txt", "a+");
    fprintf(fp,
            "input frequence =%g,phase phi=%g, Amp=%g,FFT pts n=%d, sample pts number=%d,imethod=%d,sizeof(myint)=%d,sizeof(myint2)=%u\n",
            frequence, phi,Ampltitude,n, sample_pts_num,imethod,sizeof(myint),sizeof(myint2));
    printf("FFT pts n=%d,sample pts number=%d\n", n, sample_pts_num);

   printxy(fp,n,x,y);
    printf("\n");
//如果选择复数函数计算
    if (imethod==BYFFTCOMPLEX)
    {
        for (i = 0; i < 32; i++)
        {
            f[i].real = x[i];
            f[i].imag = 0.0;
        }

    }

    switch(imethod)
    {
    case BYFFT:
        FFT(x, y, n, 0);//FFT的结果,FFT(x,y,16,0),0表示做FFT
        break;
    case BYFIXFFT:
        fix_fft(x,y,log(n)/log(2),0);//fixed point number fft
        break;
    case BYFFTfixed:
        fix_FFT(x, y, n, 0);
        break;
    case BYFFTCOMPLEX:
        fft(32,f);
        break;
    default:
        fix_fft(x,y,log(n)/log(2),0);
    }

   printxy(fp,n,x,y);
    printf("\n");

    switch(imethod)
    {
    case BYFFT:
        FFT(x, y, n,1);//IFFT的结果,FFT(x,y,16,1),1表示做IFFT
        break;
    case BYFIXFFT:
        fix_fft(x,y,log(n)/log(2),1);//fixed point number fft
        break;
    case BYFFTfixed:
        fix_FFT(x, y, n, 1);
        break;

    case BYFFTCOMPLEX:
        ifft(32,f);
        break;
    default:
        fix_fft(x,y,log(n)/log(2),1);
    }

    printxy(fp,n,x,y);

    printf("\n");

    fclose(fp);
    return 0;
}

int main()
{
    float f = 100;		//默认频率100Hz
    float phi = 1.57;	//默认初始相位pi/2,rad
    float A = 100;		//默认幅度100
    int sample_pts_num = 32;	//采样点数
    enum FFTmethod nmethod;
    short imethod=BYFIXFFT;
    int k = 1,n,m;
    int M1 = 32;
    int flag=0;
    int Amp=127;
    myint x[N],y[N];

    main2();

    printf ("input 1 to 8 corresponding  2^(k+3)=16,32,64,128,256,...,2048 points FFT,flag,Amp:\n");
    scanf("%d %d %d", &k,&flag,&Amp);

    printf("input frequence f,phase phi, imethod:\n");
    scanf("%g %g %d", &f, &phi, &imethod);
    printf("input frequence f=%g,phase phi=%g, imethod=%d,flag=%d\n", f, phi,
           imethod,flag);
    switch (k)
    {
    case 1:
    case 2:
    case 3:
    case 4:
    case 5:
    case 6:
    case 7:
    case 8:
        M1 = (1 << (k + 3));
        sample_pts_num = M1;
        break;
    case 0:
        M1 = 64;
        sample_pts_num = 32;
        break;
    case 9:
        M1 = 128;
        sample_pts_num = 64;
        break;
    default:
        M1 = 256;
        sample_pts_num = 128;
        break;
    }
    printf("M1=%d,sample_pts_num=%d\n", M1, sample_pts_num);
    if(imethod==BYPRINTSINE)
    {
        //打印8/16/24/32bit整数正弦表到文件
        for(m=1; m<5; m++)
        {
            printfixedsin(m<<3,10,flag);
        }
        for(m=7; m<12; m++)
        {
           // printfixedsin(1<<m,10,flag);
        }
    }

//循环仿真测试2**k以下点数全周期、半周期采样的FFT精度和计算时间
    for (n = 1; n<=k; n++)
    {
        fftsamples(x,y,n,1<<n,Amp,f,phi);// create samples data for fft
        if(imethod!=BYPRINTSINE)
        {
            callfft(x,y,imethod,0,n,n,Amp,f,phi);//全周期采样点
        }
    }
    for (n = 1; n<= k; n++)
    {
        fftsamples(x,y,n,1<<(n-1),Amp,f,phi);// create samples data for fft
        if(imethod!=BYPRINTSINE)
        {
            callfft(x,y,imethod,0,n,n-1,Amp,f,phi);//1/2周期采样点
        }
    }

    return 0;
}

void testsdr()
{
    //测试软件无线电
//天线或传感器
//低噪声高频放大
//ZIF零中频直接采样
//滤波:FIR或FFT
//排序或分类：前n个特征频率
//ifft
//信号输出

}

void testclassify()
{
    //测试识别分类
//天线或传感器
//低噪声高频放大
//ZIF零中频直接采样
//滤波:FIR或FFT
//排序或分类：前n个特征频率
//输出特征符号

}

void testface()
{
    //测试识别分类
//天线或传感器
//低噪声高频放大
//ZIF零中频直接采样
//滤波:FIR或FFT
//排序或分类：前n个特征频率
//输出特征符号

}

void testslam()
{
    //测试识别分类
//天线或传感器
//低噪声高频放大
//ZIF零中频直接采样
//滤波:FIR或FFT
//排序或分类：前n个特征频率
//输出定位或地图

}
